! MATRIXOP
! 
!   Implementation of several matrix-matrix products
! 
! HISTORY
! 
!   20110331 KP - Added documentation, better structure.
!   2008---- BV - Initial version
! 
! AUTHOR
! 
!   Koen Poppe, Department of Computer Science,
!   Katholieke Universiteit Leuven, Celestijnenlaan 200A,
!   B-3001 Heverlee, Belgium
!   Email:  Koen.Poppe@cs.kuleuven.be
!
!   Bart Vandewoestyne, Department of Physics,
!   Katholieke Universiteit Leuven (Kortrijk), Etienne Sabbelaan 53,
!   B-8500 Kortrijk, Belgium
!   Email:  Bart.Vandewoestyne@kuleuven-kortrijk.be
!
module matrixop
    implicit none
    private
    
    integer, public, parameter :: dp = selected_real_kind(15,307)

    ! 1. Three nested loops
    !
    ! NOTE: use the following convention for the indices
    !       i = row index of A
    !       j = column index of B
    !       k = column index of A
    !
    public :: a_maal_b_ijk
    public :: a_maal_b_ikj
    public :: a_maal_b_jik
    public :: a_maal_b_jki
    public :: a_maal_b_kij
    public :: a_maal_b_kji

    ! 2. Two nested loops with vector operations
    public :: a_maal_b_ikj_vect
    public :: a_maal_b_jki_vect
    public :: a_maal_b_kij_vect
    public :: a_maal_b_kji_vect

    ! 3. Two nested loops with dot_product
    public :: a_maal_b_ij_dot_product
    public :: a_maal_b_ji_dot_product
    
    ! 4. Two nested loops with dot_product and explicit transpose of matrix A
    public :: a_maal_b_transp_ij_dot_product
    public :: a_maal_b_transp_ji_dot_product

    ! 5. Using BLAS
    public :: a_maal_b_blas
    
    ! 6. In blocks
    public :: a_maal_b_blocks

    ! 7. Intrinsic matmul function
    public :: a_maal_b_matmul

contains

    !--------------------------------------------------------------------------
    ! 1. Three nested loops
    !--------------------------------------------------------------------------

    subroutine a_maal_b_ijk(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n   
	c = 0.0_dp ! TODO: complete this subroutine     
	n = size(A,2)
	do i=1,n
		do j=1,n
			do k=1,n
			c(i,j)=c(i,j)+a(i,k)*b(k,j)
			enddo
		enddo
	enddo
	
    end subroutine a_maal_b_ijk

    subroutine a_maal_b_ikj(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n  
        c = 0.0_dp ! TODO: complete this subroutine      
	n = size(A,2)
	do i=1,n
		do k=1,n
			do j=1,n
			c(i,j)=c(i,j)+a(i,k)*b(k,j)
			enddo
		enddo
	enddo
    end subroutine a_maal_b_ikj

    subroutine a_maal_b_jik(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n
        c = 0.0_dp ! TODO: complete this subroutine        
	n = size(A,2)
	do j=1,n
		do i=1,n
			do k=1,n
			c(i,j)=c(i,j)+a(i,k)*b(k,j)
			enddo
		enddo
	enddo
    end subroutine a_maal_b_jik

    subroutine a_maal_b_jki(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n   
        c = 0.0_dp ! TODO: complete this subroutine     
	n = size(A,2)
	do j=1,n
		do k=1,n
			do i=1,n
			c(i,j)=c(i,j)+a(i,k)*b(k,j)
			enddo
		enddo
	enddo
    end subroutine a_maal_b_jki

    subroutine a_maal_b_kij(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n     
        c = 0.0_dp ! TODO: complete this subroutine   
	n = size(A,2)
	do k=1,n
		do i=1,n
			do j=1,n
			c(i,j)=c(i,j)+a(i,k)*b(k,j)
			enddo
		enddo
	enddo
    end subroutine a_maal_b_kij
    
    subroutine a_maal_b_kji(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n    
        c = 0.0_dp ! TODO: complete this subroutine    
	n = size(A,2)
	do k=1,n
		do j=1,n
			do i=1,n
			c(i,j)=c(i,j)+a(i,k)*b(k,j)
			enddo
		enddo
	enddo
    end subroutine a_maal_b_kji

    !--------------------------------------------------------------------------
    ! 2. Two nested loops with vector operations
    !--------------------------------------------------------------------------

    subroutine a_maal_b_ikj_vect(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n   
	c = 0.0_dp ! TODO: complete this subroutine     
	n = size(A,2)
	do i=1,n
		do k=1,n
			!do k=1,n
			!c(i,j)=c(i,j)+a(i,k)*b(k,j)
			!enddo
			c(i,:) = c(i,:)+a(i,k)*b(k,:)
		enddo
	enddo
    end subroutine a_maal_b_ikj_vect

    subroutine a_maal_b_jki_vect(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n   
	c = 0.0_dp ! TODO: complete this subroutine     
	n = size(A,2)
	do j=1,n
		do k=1,n
			!do k=1,n
			!c(i,j)=c(i,j)+a(i,k)*b(k,j)
			!enddo
			c(:,j) = c(:,j)+a(:,k)*b(k,j)
		enddo
	enddo
    end subroutine a_maal_b_jki_vect

    subroutine a_maal_b_kij_vect(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n   
	c = 0.0_dp ! TODO: complete this subroutine     
	n = size(A,2)
	do k=1,n
		do i=1,n
			!do k=1,n
			!c(i,j)=c(i,j)+a(i,k)*b(k,j)
			!enddo
			c(i,:) = c(i,:)+a(i,k)*b(k,:)
		enddo
	enddo
    end subroutine a_maal_b_kij_vect

    subroutine a_maal_b_kji_vect(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n   
	c = 0.0_dp ! TODO: complete this subroutine     
	n = size(A,2)
	do k=1,n
		do j=1,n
			!do k=1,n
			!c(i,j)=c(i,j)+a(i,k)*b(k,j)
			!enddo
			c(:,j) = c(:,j)+a(:,k)*b(k,j)
		enddo
	enddo
    end subroutine a_maal_b_kji_vect

    !--------------------------------------------------------------------------
    ! 3. Two nested loops with dot_product
    !--------------------------------------------------------------------------

    subroutine a_maal_b_ij_dot_product(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n   
	c = 0.0_dp ! TODO: complete this subroutine     
	n = size(A,2)
	do i=1,n
		do j=1,n
			c(i,j) = dot_product(a(i,:),b(:,j))
		enddo
	enddo
    end subroutine a_maal_b_ij_dot_product

    subroutine a_maal_b_ji_dot_product(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer i,j,k,n   
	c = 0.0_dp ! TODO: complete this subroutine     
	n = size(A,2)
	do j=1,n
		do i=1,n
			c(i,j) = dot_product(a(i,:),b(:,j))
		enddo
	enddo
    end subroutine a_maal_b_ji_dot_product

    !--------------------------------------------------------------------------
    ! 4. Two nested loops with dot_product and explicit transpose of matrix A
    !--------------------------------------------------------------------------

    subroutine a_maal_b_transp_ij_dot_product(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	real(kind=dp), dimension(size(a,2),size(a,1)) :: at
	integer i,j,k,n	
	c = 0.0_dp ! TODO: complete this subroutine  
	at = transpose(a)   
	n = size(A,2)
	do i=1,n
		do j=1,n
			c(i,j) = dot_product(at(:,i),b(:,j))
		enddo
	enddo
    end subroutine a_maal_b_transp_ij_dot_product

    subroutine a_maal_b_transp_ji_dot_product(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	real(kind=dp), dimension(size(a,2),size(a,1)) :: at
	integer i,j,k,n	
	c = 0.0_dp ! TODO: complete this subroutine  
	at = transpose(a)   
	n = size(A,2)
	do j=1,n
		do i=1,n
			c(i,j) = dot_product(at(:,i),b(:,j))
		enddo
	enddo
    end subroutine a_maal_b_transp_ji_dot_product

    !--------------------------------------------------------------------------
    ! 5. Using BLAS
    !--------------------------------------------------------------------------

    subroutine a_maal_b_blas(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
	integer n
        c = 0.0_dp ! TODO: complete this subroutine
	n = size(A,2)
	call dgemm('N','N',n,n,n,1.0_dp,a,n,b,n,0.0_dp,c,n)
    end subroutine a_maal_b_blas
    
    !--------------------------------------------------------------------------
    ! 6. In blocks
    !--------------------------------------------------------------------------

    subroutine a_maal_b_blocks(a, b, c, d)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
        integer, intent(in) :: d
	integer i,j,k,n, blocks
	n= size(A,2)	
	blocks = n/d
	c = 0.0_dp ! TODO: complete this subroutine     
	do i=0,blocks-1
		do j=0,blocks-1
			do k=0,blocks-1
			c((i*d+1):(i+1)*d,(j*d+1):(j+1)*d)=c((i*d+1):(i+1)*d,(j*d+1):(j+1)*d)+&
			matmul(a((i*d+1):(i+1)*d,(k*d+1):(k+1)*d),b((k*d+1):(k+1)*d,(j*d+1):(j+1)*d))
			enddo
		enddo
	enddo
	
    end subroutine a_maal_b_blocks

    !--------------------------------------------------------------------------
    ! 7. Intrinsic matmul function
    !--------------------------------------------------------------------------

    subroutine a_maal_b_matmul(a, b, c)
        real(kind=dp), dimension(:,:), intent(in)  :: a, b
        real(kind=dp), dimension(:,:), intent(out) :: c
        c = matmul( a, b ) ! Already completed
    end subroutine a_maal_b_matmul

end module matrixop
